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Abstract 

We study static, spherically symmetric solutions in a recently proposed ghost-free model of non-linear 
massive gravity. We focus on a branch of solutions where the helicity-0 mode can be strongly coupled 
within certain radial regions, giving rise to the Vainshtein effect. We truncate the analysis to scales 
below the gravitational Compton wavelength, and determine analytically the number and properties 
of local solutions which exist asymptotically on large scales, and of local (inner) solutions which exist 
on small scales. We find two kinds of asymptotic solutions, one of which is asymptotically flat, while 
the other one is not. We find also two types of inner solutions, one of which displays the Vainshtein 
mechanism, while the other exhibits a self-shielding behaviour of the gravitational field. We analyse in 
detail in which cases the solutions match in an intermediate region. The asymptotically flat solutions 
connect only to inner configurations displaying the Vainshtein mechanism, while the non asymptotically 
flat solutions can connect with both kinds of inner solutions. We show furthermore that there are some 
regions in the parameter space where global solutions do not exist, and characterise precisely in which 
regions of the phase space the Vainshtein mechanism takes place. 
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1 Introduction 



Promoting Einstein's gravity to a classical theory of a massive graviton is a theoretical challenge. The 
appearance of a ghost instability and the fact that in the massless limit solutions may not reduce to those 
of General Relativity (GR) have been the two main obstacles for a successful Lagrangian construction [1]. 
Both of these problems are related to the scalar sector of the theory. In recent years, new developments 
have opened a window into a ghost-free model for massive gravity [2, 3, 4], based on the Galilean symmetry 
for the scalar mode [5]. However, with respect to the second issue the story is not entirely settled, as we 
will discuss here. 

In this massive gravity theory, spherically symmetric solutions do not obey a uniqueness principle, 
like Birkhoff's theorem in General Relativity. As a result, the theory exhibits two classes of static spher- 
ically symmetric solutions, which differ in how the scalar sector couples to graviton. The first class of 
solutions naturally contains (Anti) de Sitter-Schwarzschild solutions, where the cosmological constant is 
proportional to the graviton's mass [6, 7, 8]. Some of these solutions are very similar to those found 
many years ago in the simple Fierz-Pauli model [9] , and what characterises them is that the scalar mode 
is strongly coupled everywhere, which avoids a discontinuous connection to GR in the massless limit. 
Therefore, at the level of the background, these solutions are indistinguishable from their counterparts 
in General Relativity, but do present differences at first orders in perturbations. Interestingly, in some 
of these solutions, scalar and vector perturbations are not dynamical at first order in perturbations, and 
only the tensor modes differ significantly from General Relativity [10, 11, 12]. This suggests that only 
the detection of gravitational waves can rule out these solutions in favour of ACDM. 

In the second class of solutions the story is not that clear. The scalar sector is not strongly coupled 
everywhere, and the solutions behave very differently as a function of the radial coordinate. The purpose 
of this paper is to exhibit the rich structure of such class of solutions, and show whether one can recover 
GR solutions in the massless limit, at least within some radial region. This phenomenon of recovering GR 
predictions within some macroscopic radius from a mass source is known as the Vainshtein mechanism 
[13]. In the range where the solution reduces to GR the scalar field becomes strongly coupled. Therefore, 
it is imperative to include non-linearities in the scalar sector in order to study the Vainshtein mechanism. 
Solutions presenting these behaviour in massive gravity where found in [14], and more recently in the 
ghost-free model considered here in [6, 7, 15, 16]. Unfortunately, these solutions do not cover the whole 
parameter space of the theory as we will describe here. 

We organise the paper as follows: in Section 2 we give an overview of the main equations describing the 
theory and spherically symmetric solutions, with particular attention to this second branch. In Section 3 
we present an exhaustive analysis of the vacuum equations, with particular attention to the behaviour of 
their solutions near the origin and towards infinity. In Section 4 we present our main results, including 
the parameter space in which GR solutions are recovered via the Vainshtein mechanism and numerical 
solutions of the vacuum equations for some representative cases. Finally, we present some conclusions in 
Section 5. 

2 Ghost-free massive gravity and spherically symmetric ansatze 

We consider the following Lagrangian for massive gravity, a non-linear extension of Fierz-Pauli theory 
proposed in [3] 




(2.1) 
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The potential U depends on a dimension-full parameter m, which sets the graviton mass scale, and on 
two dimensionless parameters 03 and a 4 . It has the following functional form 

U = -m 2 [U 2 + a 3 U 3 + a 4 U 4 ] , (2.2) 

with 

U 2 = (tr/C) 2 -tr(/C 2 ), 

U 3 = (tr/C) 3 - 3(tr/C)(tr/C 2 ) + 2tr/C 3 , 

U A = (tr/C) 4 - 6(tr/C) 2 (tr/C 2 ) + 8 (tr/C) (tr/C 3 ) + 3(tr/C 2 ) 2 - 6tr/C 4 . 
The tensor is defined as [3] 

K» = ^-(^S)^ , (2.3) 

where the square root of a tensor is formally understood as \fM. ?VM-a = M.^, for any tensor . 
The tensor T,^ is a fiducial metric, which also makes the theory reparametrisation invariant by means of 
four scalars </> M , and it is given by 

5V = d^ a d„4> b Vab . (2.4) 

The theory defined by (2.1) has Minkowski spacetime as solution, hence one can rewrite the metric g^ u 
and the scalars ^ as deviations from flat space, namely 

g^u = r)n» + V. ^ = xli + ^ ( 2 - 5 ) 

where x M are the usual cartesian coordinates spanning r]^ u . Therefore, a change of coordinates — > x^+^ 
should be accompanied by the following transformation of the Stiickelberg field 7r M , 

^ -»• vr^ + £ M , (2.6) 

in order to recover full diffeomorphism invariance. Moreover, we choose the unitary gauge, where 7r M = 0, 
and the potential (2.2) considerably simplifies. 

Since we are intested in spherically symmetric solutions, we use the most general static ansatze, given 

by 

ds 2 = -C(r) dt 2 + A{r) dr 2 + 2D(r) dtdr + B(r)dn 2 , (2.7) 

where d£l 2 = d# 2 +sin 2 Odtfi 2 . We choose to write the non-dynamical flat metric as ds 2 = —dt 2 +dr 2 +r 2 dQ 2 . 
It should be noticed that this is not a coordinate choice, but a way to simplify the expressions. We plug 
the previous metric into the Einstein equations = Tjf v , where the energy momentum tensor is defined 

as = -^L 5 ^~® V U . The Einstein tensor satisfies the identity D(r) G tt + C(r) G tr = 0, which 

implies the algebraic constraint = D(r) T^f + C(r) Tj£. This last equation is solved in two possible ways, 
defining two class of solutions: either the metric is diagonal D = 0, or B = b^r 2 , where 60 is a function 
of Q3 and 04 only. Notice that this classification only holds in the unitary gauge, since one can always 
map the metric from one class to the other by a coordinate transformation, but to the price of exciting 
components of n^. 

The class of solutions with a non-diagonal metric leads to Schwarzschild or Schwarzschild-de Sitter 
solutions [6, 7, 8], as explained in the introduction. In this sector of the theory GR is recovered by means 
of a strongly coupled 7T M everywhere [10, 11]. However, in the other class of solutions, where the metric is 
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diagonal in the unitary gauge, the situation is different. As we discuss in the next section, 7r^ may or may 
not be strongly coupled: it could be strongly coupled within certain radial region, leading to a Vainshtein 
effect. This branch is the one that concerns us here, so from now on we will only consider static diagonal 
spherically symmetric solutions in the unitary gauge. 

The problem of finding exact vacuum solutions in this branch is an open question, but one can make 
progresses by considering perturbations (not necessarily small) from flat space, and the following ansatz 
results adequate for this purpose, 



ds 2 



(l + iV(r)) 2 dt 2 + (l + F(r)) V 2 + r 2 {l + H{r)^ ' dfi 2 . (2.8) 



In order to analyse the system, it is convenient to introduce a new radial coordinate 

" = TTWy (2 ' 9) 

so that the linearised metric is expressed as 

ds 2 = -(1 + n)dt 2 + (1 - f)d P 2 + p 2 dn 2 , (2.10) 

where f(p) = F(r(p)) - 2h(p) - 2pti(p), n(p) = 2N(r(p)), h(p) = H{r(p)) and a prime denotes a 
derivative with respect to p. As discussed above, one should be careful with this change of coordinates, 
since, after fixing a gauge, a change of frame in the metric modifies the Stiickelberg field tt p as well. 
It turns out that this coordinate transformation excites the radial component of n^, which explicitly is 
ir p = ph. Therefore, from now on one can think of h as simply being the only non-zero component of 
the Stiickelberg field ir^-. At linear order, the equations for the functions n(p), f(p) and h(p) in the new 
variable p are 

= (m 2 p 2 + 2) f + 2p(f + m 2 p 2 ti + 3m 2 ph) , (2.11) 
= -m 2 p 2 (n-4h) - pri - f, (2.12) 

= f+ l -pn'. (2.13) 

In this linear expansion, the solutions for n and / are 

„ _ J™—, 

3p 
AGM 

f = ~(l + m/))e-T (2.14) 

3/9 

where we fix the integration constant so that M is the mass of a point particle at the origin, and 8ttG = 
M~ 2 . These solutions exhibit the vDVZ discontinuity, since the post-Newtonian parameter 7 = f/n is 
7 = i(l + mp), which in the massless limit reduces to 7 = 1/2, in disagreement with GR and Solar system 
observations. 

However, in order to understand what really happens in this limit, we must also analyse the behaviour 
of h, or equivalently ir p , as m — > 0. To do this, we consider scales below the Compton wavelength rap <C 1, 
and at the same time ignore higher order terms in GM. Under these approximations, the equations of 
motion can still be truncated to linear order in / and n, but since h is not necessarily small, we have to 
keep all non-linear terms in h. In other words, we take the usual weak field limit for the metric fields, 
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but keep all non-linearities in the Stiickelberg field, since we expect regions where this field is strongly 
coupled. As shown in [7], the field equations reduce to the following system of coupled equations for the 
fields /, n, h: 

f = -2^y- ~ {mpf [h-(l + 3a 3 )h 2 + (a 3 + 4q 4 )/i 3 ] (2.15) 



1CM 

pn' = — p (mp) 2 [h - (a 3 + ia 4 )h 3 } (2.16) 

— [1 - 3(a 3 + 4a 4 )/i 2 ] = -{mpfi^h - 3(1 + 3a 3 )h 2 + [(1 + 3q 3 ) 2 + 2(q 3 + 4a 4 )] h 3 

P 3 (2 ' 17) 

- -(a 3 + 4a 4 ) 2 /i 5 j 



These equations can also be obtained directly from the decoupling theory [2], as it was shown in [7]. The 
previous expressions are the starting point of our analysis and classifications of solutions, which we will 
discuss in the next Section. 



3 Classifying Solutions 

Our aim is to scan the (a 3 ,a 4 ) parameter phase space of theories, to understand how many solutions the 
system (2.15)-(2.17) admits, and characterise their (asymptotic) geometrical properties. 

Since the last equation (2.17) does not contain / and n, the field h obeys a decoupled equation, thus 
once a solution of this equation is found, the gravitational potentials /, n are uniquely determined (up 
to an integration constant) by the other two equations (2.15) and (2.16). Therefore, we first focus on 
classifying the number and properties of solutions of (2.17) in every point of the phase space, and then 
discuss the behaviour of the gravitational potentials correspondent to these solutions. 

3.1 The quintic equation 

For notational convenience, we define a = 1 + 3 a 3 and f3 = a 3 + 4 a 4 . Therefore, the system that 
determines the gravitational potentials in terms of h takes the form 

f = _2^--(mp) 2 (h-ah 2 + ph 3 ) (3.1) 

n > = 2 ( ^--m 2 p(h-ph 3 ) (3.2) 
while the equation which determines h reduces to 

| P 2 h 5 (p) - (a 2 + 2/3) h 3 {p) + 3 (a + (3A(p)) h 2 (p) - | h{p) - A(p) = (3.3) 

where A(p) = (p v /p) 3 and p v is the Vainshtein radius defined as p v = (GM/m 2 ) 1 ^ 3 . The last equation 
is an algebraic equation in h, A, a and (3; at fixed p, a and (3 it is, in fact, a polynomial equation of 
fifth degree in h, except for the special case f3 = 0. In this particular case of /3 = 0, the equation for h 
becomes a cubic equation and it is possible to obtain solutions for h and the metric perturbations exactly. 
These solutions were studied in [7] and it was shown that the solutions exhibit the Vainshtein mechanism. 
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Therefore, in what follows, we assume /3 7^ 0. It is difficult to find analytical solutions to this quintic 
equation, and it is not easy even to understand how many solutions it admits. In fact, even if a local 
solution is found in a radial interval, it is not, in general, possible to extend it to the whole radial domain, 
as we will explicitly see in the next section. 

Nevertheless, it is possible to determine exactly how many local solutions exist in a neighbourhood 
of p = +cx), which we refer as asymptotic solutions, and also how many local solutions exist in a neigh- 
bourhood of p = + , which we call inner solutions. Furthermore, we can find analytically their leading 
behaviour as a function of p. Any global solution of (3.3) should necessarily interpolate between one of 
the asymptotic solutions and one of the inner solutions. Therefore, our aim is to understand, for each 
point in the (a, (3) phase space, whether and how the above solutions match. 

A systematic approach to vainshtein effects in covariant galileon theory was performed in [17], and 
in general scalar-tensor theories [18]. However, in massive gravity a systematic treatment was not fully 
performed, missing some asymptotic behaviours and assuming the matching between inner solutions and 
the asymptotic ones always occurs. 

3.2 Phase space analysis 

To be able to describe how the matching works in all the phase space, in principle we should study 
separately every point (a, (3). However, this is not necessary since equation (3.3) obeys a remarkable 
symmetry: defining the quintic function as 

q (h; A, a, 0) = | /3 2 h 5 - (a 2 + 2/3) h 3 + 3 (a + pA) h 2 - - h - A (3.4) 
it is simple to see that 

q^;^,ka,k 2 (3) =^q(h;A,a,l3) (3.5) 

Therefore if a local solution of (3.3) exists for a given (a, /3) within a certain radial interval, it would also 
be present for (ka, k 2 /3), for k > 0, with h being replaced by h/k and the radial interval rescaled by 1/^fk. 
As a result, each point belonging to the a > part of the parabola /3 = ca 2 (with c any non- vanishing 
constant) of the phase space shares the same physics, hence having the same number of global solutions 
and matching properties. The same is true for the points belonging the a < part of the parabola. So, 
to understand the global structure of the phase space, it is sufficient to analyse one point for each of the 
half parabolas present in the phase space. 

It is worthwhile to point out that our starting equations (3.1)-(3.3) were constructed assuming GM < 
p < 1/m, but in the following analysis we use the whole radial domain < p < +00 . On one hand, this 
allows us to characterise exactly the number and properties of solutions on large and small scales. On 
the other hand, the picture we have in mind is that the Compton wavelength of the gravitational field 
p c = 1/m is of the same order of the Hubble radius today, and that there is a huge hierarchy between 
p c and the gravitational radius 1 p g = GM, i.e. p c j p g ^> 1. Therefore, we expect that extending the 
analysis to the whole radial domain captures the correct physical results. 

3.3 Asymptotic and inner solutions 

We sum up here the results obtained in the appendix on the existence and properties of asymptotic and 
inner solutions of eq. (3.3). We refer the reader to the appendix for the details. 

1 We are using units where the speed of light speed has unitary value. 
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3.3.1 Asymptotic solutions 

In a neighbourhood of p — > +00 there are, depending on the value of (a,/3), three or five solutions to 
eq. (3.3). In particular: 

- There is always a decaying solution, which we indicate with L. Its asymptotic behaviour is 

Hp) = '+ R 0>) ( 3 - 6 ) 

where linip^+oo p 3 R(p) = 0. This solution corresponds to a spacetime which is asymptotically flat, 
as one can see from eqs. (3.1)-(3.2). 

- Additionally, there are two or four solutions to eq. (3.3) which tend to a finite, nonzero value as 
p —7- +00. We name these solutions with C+, C_, Pi and P2 (details about this denomination are 
given in the appendix). Their asymptotic behaviour is 

h{p) = C+R{p) (3.7) 

where \\m p ^ +00 R(p) = and C is a root of the reduced asymptotic equation (A. 5). From eqs. (3.1)- 
(3.2), one can get convinced that these solutions correspond to spacetimes which are asymptotically 
non-flat. Interestingly, the leading term in the gravitational potentials scales as p 2 for large radii, 
the same scaling which we find in a de Sitter spacetime. 



3.3.2 Inner solutions 

In a neighbourhood of p — > + there are either one or three solutions to eq. (3.3). For f3 > there are 
exactly three inner solutions, while for j3 < there is only one inner solution. In particular: 

- There is always a diverging solution, which we denote by D. Its leading behaviour is 

KP) = ~ {f\- + R{P) (3-8) 

V p p 

where lim /9 _ > o+ (R(p) I p) i s finite. This solution exists for both (3 > and /3 < 0, with opposite signs 
for each case. Using this solution in eqs. (3.1)-(3.2), one realises that the /i 3 term cancels the GM/p 
term, so the gravitational field is self-shielded and does not diverge as p — > + . This solution is in 
strong disagreement with gravitational observations. 

- For (3 > 0, there are two additional solutions to eq. (3.3), which tend to a finite, non-zero value as 
p — > + . We indicate these solutions by F + and F_ . Their leading behaviour is 

Hp) = ±^ + R{p) (3.9) 

where lim p _ >0 + R = 0. Notice that for (3 < there are no solutions to eq. (3.3) which tend to a 
finite value as p — > + . 

The expressions (3.1)-(3.2) for the gravitational potentials imply that the metric associated to these 
solutions (F + and F_) approximate the linearised Schwarzschild metric as p — > + . 

From the behaviour of the inner solutions, one concludes that only in the f3 > part of the phase space 
solutions may exhibit the Vainshtein mechanism, but not necessarily for all values of a. In the next 
subsection we see more in detail how this mechanism works. 



7 



3.3.3 Vainshtein mechanism 



In order to study where in the phase space the Vainshtein mechanism works, it is useful to compare the 
gravitational potentials / and n with their counterparts in the GR case. In the weak field limit, the 
Schwarzschild solution of GR reads 

ds i = - (1 - ^) dt 2 + (1 + ^j^j dp 2 + p 2 dn 2 (3.10) 

so by calling /gr = uqr = —2GM/p we obtain 

/ 



Igr 

n' 
n GR 



^U^V"" 3 ) (3 - i2) 



Let us now first discuss the asymptotic solutions. For the decaying solution L, we have that the linear 
contribution in h rescales the coefficients of the Schwarzschild- like terms, so we obtain f / /gr — > 2/3 
and n'/n GR — > 4/3 for p — > +oo. For the non-decaying solutions C± and Pi,2 ; the leading behaviour 
for f / fcR an d n '/ n GR ^ s proportional to (p/pv) 3 in both cases, however the proportionality coefficients 
generally differ since they have a different functional dependence on a and f3. There are some special 
cases for (a,/3) where these asymptotic solutions lead to f/n — > 1 as p — > +oo, and therefore have the 
same behaviour as in a de Sitter spacetime. 

Consider instead the inner solutions. For the finite solutions F± we obtain (/ / /gr) — > 1 and 
(n' /n' GR ) — > 1 as p — > + , where the corrections scale like p 3 . On the contrary, for the diverging so- 
lution D, the cubic terms in h cancel out the contribution coming from the Schwarzschild- like terms, as 
explained above, and so (///gr) — > and (n'/n GR ) — > when p — > + . In this case, corrections are 
linear in p. 

Therefore, any global solution of equation (3.3) which interpolates between L and F± provides a 
realisation of the Vainshtein mechanism in an asymptotically flat spacetime, whereas an interpolation 
between C± or Pi^ with F± exhibits the Vainshtein mechanism in an asymptotically non-flat spacetime. 
Furthermore, notice that any asymptotic solution which interpolates with the inner solution D does no 
lead to the Vainshtein mechanism. These matchings will be explicitly exposed in the next section. 



4 Results 

4.1 Solutions matching 

The phase space diagram which displays our results about solution matching is given in figure 1. We 
discuss separately the /3 > and (3 < part of the phase space, and refer to the figure for the numbering 
of the regions. The notation I -H- A means that there is matching between the inner solution I and the 
asymptotic solution A. 

P < 

In this part of the phase space, there is only one inner solution, D, so there can be at most one global 
solution to (3.3). There are three distinct regions which differ in the way the matching works: 
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Figure 1: Phase space diagram in (a,/3) for the solutions to the quintic equation (3.3) in h, where the 
different regions show different matching of inner solutions to asymptotic ones. The lines splitting the 
regions are half parabolas (/3 oc a 2 , with a > or a < 0) due to rescaling symmetry of eq. (3.3). 
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- region 1: D -H- C + . In this region, there are three or five asymptotic solutions, and only one of 
them, C+, is positive. This solution is the one which connects with the inner solution D, which is 
also positive, leading to the only global solution of eq. (3.3). The boundaries of this region are the 
line f3 = for a < and the parabola (3 = cyi a 2 for a > 0, where cyi is the negative 2 root of the 
equation -4 - 8y + 88 y 2 - 1076 y 3 + 2883 y 4 = (approximatively, c 12 ~ -0.1124). 

- region 2: No matching. In this region there are three asymptotic solutions. However, none of 
them can be extended all the way to p — > + , and so, despite the fact that local solutions exist both 
at infinity and near the origin, equation (3.3) does not admit any global solution. The boundaries of 
this region are the parabola (3 = c\2 a 2 and the (negative) five-roots-at-infinity parabola (3 = c_ a 2 , 
where c_ is the only real root of the equation 8 + 48 y — 435 y 2 + 676 y 3 = (approximatively, 
c_ ~ -0.0876). 

- region 3: D <H> P2. This region coincides with the a > 0, (3 < part of the five roots at infinity 
region of the phase space (see fig. 7). The largest positive asymptotic solution, P2, is the one which 
connects to D, leading to the only global solution of eq. (3.3). 

p > 

In this part of the phase space, there are three inner solutions, D, F + and F_, so there can be at 
most three global solutions to eq. (3.3). There are six distinct regions with different matching properties: 

- region 4: F_ <H> L , D <H> C_. This region lies inside the a > 0, j3 > part of the five roots at 
infinity region of the phase space (see fig. 7), so there are five asymptotic solutions. Of the five 
asymptotic solution, C_ and L can always be extended to p — > + , while C+, Pi and P2 cannot. 
So there are just two global solutions to eq. (3.3). The boundaries of this region are the parabola 
f3 = c 45 a 2 , where c 45 = 1/12 ~ 0.0833, and the line (3 = 0. 

- region 5: F + <H> C+ , F_ <H> L, D <H> C . In this region there are three or five asymptotic solutions; 
C , C + and L can always be extended to p — > + , while Pi and P2 , where present, cannot. So 
there are three global solutions to (3.3). The boundaries of this region are the parabola (3 = C45 a 2 
for a > and the parabola (3 = c$q a 2 for a < 0, where C56 = (5 + VT3)/24 ~ 0.3586. 

- region 6: D -H- C_ , F + -H- C + . In this region there are three asymptotic solutions, however only two 
of them can be extended to p — > + , while L cannot. Therefore, there are just two global solutions 
to eq. (3.3). The boundaries of this region are the parabolas (3 = c^a 2 and (3 = cq-j a 2 , where 
C67 is the positive root of the equation — 4 — 8y + 88 y 2 — 1076 y 3 + 2883 y 4 = (approximatively, 
c 67 ~ 0.3423). 

- region 7: F + <H> C + . In this region there are three asymptotic solutions, however only one of them 
can be extended to p — > + , while L and C_ cannot. The boundaries of this region are the parabola 
(3 = C67 a 2 and the (positive) five-roots-at-infinity parabola (3 = c + a 2 , where c + = 1/4. Note that 
on the (q < 0) part of the parabola (3 = 1/3 a 2 there is the additional matching F_ o C_, so for 
these points there are two global solutions to eq. (3.3). 

- region 8: F + C + , F_ o Pj , D H P2. This region lies inside the a < 0, (3 > part of the five 
roots at infinity region of the phase space (see fig. 7), so there are five asymptotic solutions. Only 

2 The equation — 4 — 8 y + 88 y 2 — 1076 j/' 5 + 2883 y 4 — has only two real roots, one positive and one negative. 
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three of them can be extended to p — > + , while C_ and L cannot. The boundaries of this region 
are the parabolas ft = c + a 2 and ft = Cs9 a 2 , where egg = (5 — VT3)/24 ~ 0.0581. 

- region 9: F_ o Pj , D o P2. This region lies inside the a < 0, ft > part of the five roots 
at infinity region of the phase space (see fig. 7), so there are again five asymptotic solutions. The 
matching is similar to that of region 8, apart from the fact that C+ cannot be extended to p — > + 
anymore; hence there are just two global solutions to eq. (3.3). The boundaries of this region are 
the parabola ft = Cs9 a 2 and line ft = 0. 

We note that the decaying solution L never connects to the diverging one D, so we cannot have a spacetime 
which is asymptotically flat and exhibit the self-shielding of the gravitational field at the origin. On the 
other hand, finite non-zero asymptotic solutions (C± or Pi^) can connect to both finite and diverging 
inner solutions. Therefore, one can have an asymptotically non-flat spacetime which presents self-shielding 
at the origin, or an asymptotically non-flat spacetime which tends to Schwarzschild spacetime for small 
radii. More precisely, for ft < there are only solutions displaying the self-shielding of the gravitational 
field, apart from region 2 where there are no global solutions. Therefore the Vainshtein mechanism never 
works for ft < 0. In contrast, for ft > all three kinds of global solutions are present. Solutions with 
asymptotic flatness and the Vainshtein mechanism are present in regions 4 and 5, while solutions which 
are asymptotically non-flat and exhibit the Vainshtein mechanism do exist in all (ft > 0) regions but 
region 4. Finally, solutions which display the self-shielding of the gravitational field are present in all 
(ft > 0) regions but region 7. 

4.2 Numerical solutions 

We present here the numerical solutions for the h field and the gravitational potentials in some represen- 
tative cases. We choose a specific realisation for each of the three physically distinct cases, namely asymp- 
totic flatness with Vainshtein mechanism, asymptotically non-flat spacetime with Vainshtein mechanism, 
and asymptotically non-flat spacetime with self-shielded gravitational field at the origin. In addition, we 
consider the case in which there are no global solutions to eq. (3.3). This provides an illustration of what 
happens, in general, to local solutions of eq. (3.3) which cannot be extended to the whole radial domain, 
and give an insight on the phenomenology of equation (3.3). 

Asymptotic flatness with Vainshtein mechanism 

Let's consider the case in which the solution of eq. (3.3) connects to the decaying solution at infinity 
L and to a finite inner solution (in this case F_). In figure 2, the numerical solutions for h (dashed 
line), / / fcR (bottom continuous line) and n' /n' GR (top continuous line) are plotted as functions of the 
dimensionless radial coordinate x = p/p v . These solutions correspond to the point (a, ft) = (0,0.1) of 
the phase space. 

This plot displays very clearly the presence of the vDVZ discontinuity and its resolution via the 
Vainshtein mechanism. For large scales, h is small and the gravitational potentials behave like the 
Schwarzschild one, however their ratio is different from one, unlike the massless case. Note that the 
ratio of the two potentials for p 3> p v is independent of m, so does not approach one as m — > (vDVZ 
discontinuity). However, on small scales h is strongly coupled, and well inside the Vainshtein radius 
the two potentials scale again as the Schwarzschild one, but their ratio is now one even if m 7^ 0. So, 
the strong coupling of the h field on small scales restores the agreement with GR (Vainshtein mechanism) . 
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Figure 2: Numerical solutions for the case F_ f-> L. 



Asymptotically non-flat spacetime with Vainshtein mechanism 

Let's consider now the case in which the solution of eq. (3.3) connects to a finite solution at infinity 
and to a finite inner solution. We consider for definiteness the phase space point (a,/3) = (0,0.1). In 
figure 3, we plot the numerical results for the gravitational potentials (normalised to their GR values) 
and the global solution of eq. (3.3) which interpolates between the inner solution F + and the asymptotic 
solution C+. 




0.0 0.5 1.0 1.5 2.0 2-5 



Figure 3: Numerical solutions for the case F + C + . 
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We can see that, on large scales, the gravitational potentials are not only different one from the other 
but also behaves very differently compared to the GR case. However, on small scales there is a macro- 
scopic region where the two potentials agree, and their ratio with the Schwarzschild potential stays nearly 
constant and equal to one. Therefore, also in this case the small scale behaviour of h guarantees that GR 
results are recovered, even if the spacetime is not asymptotically flat. This behaviour provide then, in a 
more general sense, a realisation of the Vainshtein mechanism. 

Asymptotically non-flat spacetime with self-shielding 

We turn now to the case where the solution of eq. (3.3) connects to a finite solution at infinity and 
to the diverging inner solution. In figure 4, we plot the global solution h and the associated gravitational 
potentials, normalised to their GR values, correspondent to the phase space point (ct,f3) = ( — 1,-0.5). 
It is apparent that there are no regions where the solutions behave like in the GR case. 




Figure 4: Numerical solutions for the case D f-> C+. 



To see that the gravitational potentials are indeed finite at the origin, we plot in figure 5 the potentials 
/ and n' themselves, as functions of p/p v . We choose for definiteness the following ratio between the Comp- 
ton wavelength and the gravitational radius p c j ' p g = 10 6 , and plot the potentials for 0.01 < p/p v < 2. 
Note that, since in this case p c / ' p v = ypc/pg = 10 2 , the range where the functions are plotted is well 
inside the range of validity of our approximations. We can see that the potentials approach a finite value 
as p — > + , and so indeed the gravitational field does not diverge at the origin. 

No matching 

Finally, we consider the case in which equations (3.1) — (3.3) do not admit global solutions. We consider 
for definiteness the phase space point (a,/3) = (1 , —0.092). In figure 6 we plot all the local solutions of 
the quintic equation (3.3) as functions of the dimensionless radial coordinate x = p/p v . 

For < x < 0.4, there is only one local solution (the top continuous curve), which connects to the 
diverging inner solution D. At x ~ 0.4 a pair of solutions is created (dashed and continuous negative 
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Figure 6: Numerical results for all local solutions of eq. (3.3) in the case where there is no matching. 
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valued curves), and at x ~ 0.9, another pair of solutions is created (positive valued dashed curve and 
positive valued bottom continuous curve). However, at x ~ 1.3 one of the newly created functions (the 
positive valued dashed curve) annihilates with the solution with connects to the inner solution, so for 
x > 1.3 there are three local solutions, which finally connect with the asymptotic solutions C_ , L and 
C+. Therefore, the number of existing local solutions is one for < x < 0.4, three for 0.4 < x < 0.9, five 
for 0.9 < x < 1.3 and three for x > 1.3. We can see that, despite the fact that for every p there is at 
least one local solution, there does not exist a solution which extends over the whole radial domain. Note 
that the solutions are created and annihilated in pairs; furthermore, the pairs of solutions have infinite 
slope when they are created and when they annihilate, which results in the gravitational potentials having 
diverging derivatives when the creation/annihilation points are approached. These are found to be general 
features of the phenomenology of equation (3.3). In fact, in most part of the phase space some of the 
asymptotic/inner solutions cannot be extended to the whole radial domain < p < +oo. However, it 
never happens that the value of the local solutions diverges at a finite radius, while they always disappear 
or are created in pairs, with their values remaining bounded but with derivatives which diverge. 

5 Conclusions 

Recently, a ghost-free model of non-linear massive gravity has been proposed. We studied static, spheri- 
cally symmetric solutions of this theory inside the Compton radius of the gravitational field, and considered 
the weak field limit for the gravitational potentials, while keeping all non-linearities of the helicity-0 mode. 
For every point of the two free parameter phase space, we characterised completely the number and prop- 
erties of asymptotic solutions on large scales and also of inner solutions on small scales. In particular, 
there are two kinds of asymptotic solutions, where one of them is asymptotically flat and the other one 
is not. There are also two kinds of inner solutions, one which displays the Vainshtein mechanism and the 
other which exhibits the self-shielding of the gravitational field near the origin. 

We described under which circumstances the theory admits global solutions interpolating between the 
asymptotic and inner solutions, and found that the asymptotically flat solution connects only to inner 
solutions displaying the Vainshtein mechanism, while solutions which diverge asymptotically can connect 
to both kinds of inner solutions. Furthermore, we showed that there are some regions in the parameter 
space where global solutions do not exist, and characterised precisely in which regions of the phase space 
the Vainshtein mechanism is working. 

Our study embraces all of the phase space spanned by the two parameters of the theory. Notably, we 
found that, within our approximations, the asymptotic and inner solutions cannot in general be extended 
to the whole radial domain. In particular, we exhibited extreme cases in which global solutions do not 
exist at all. This happens because at a finite radius the derivatives of the metric components diverge, 
while the metric components themselves remain bounded. When the derivatives of the metric cease to 
be small, the approximations we used to derive the equations (3.1) — (3.3) break down. It would be 
interesting to study what happens at this radius in the full theory. 
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A Asymptotic behaviour of solutions of the quintic equation 



In this appendix we study the local solutions of the quintic equation (3.3) in a neighbourhood of p = + 
and also in a neighbourhood of p = +00. As previously mentioned we consider only the /3 7^ case, since 
for (3 = the equation becomes a cubic and it is possible to find all local solutions analytically (see [7]). 

A.l Asymptotic and inner solutions 
Asymptotic solutions 

Suppose that a solution h(p) of the quintic equation (3.3), which for convenience we rewrite here, 

I f3 2 h 5 - (a 2 + 2/3) h 3 + 3 (a + f3A^j h 2 - | h - A = (A.l) 

exists in a neighbourhood of p = +00 , and that it has a well defined limit as p — > +00. Then such a 
solution cannot be divergent. To see this, suppose that indeed the solution is divergent | linip^+oo h(p)\ = 
+00 : it is then possible to divide the quintic equation by h 5 (in a neighbourhood of p = +00 ), obtaining 
the following equation in v = 1/h 

If _ ( a 2 + 2/3 ) w 2 + 3 ( a + 0^) „3 _ 3 y 4 _ Av 5 = (A 2) 

Taking the p — >■ +00 limit of both sides of this expression one obtains j3 = 0, which is precisely against 
our initial assumption. 

Suppose now that lim.p^ +00 h(p) is finite, and let's call it C. Then both of the sides of the quintic 
equation itself have a finite limit when p — > +00 , and taking this limit one gets 

^/3 2 C 5 - (a 2 + 2/3) C 3 + 3aC 2 -^C = (A.3) 

It follows then that the allowed asymptotic values at infinity for h(p) are the roots of the following 
equation, which we call the asymptotic equation 

s/(y) = I /3 2 y 5 - {a 2 + 2/3) y 3 + 3 a y 2 - | y = (A.4) 

Note that y = is always a root of this equation, and in fact a simple root (i.e. a root of multiplicity one) 
since Jj^(O) ^ . Dividing by y, one obtains that the other asymptotic values for h(p) are the roots of 
the reduced asymptotic equation 

<(y) = ^(3 2 y 4 - (a 2 + 2p)y 2 + 3ay-^ = (A.5) 

This last equation is a quartic, so it can have up to 4 (real) roots, depending on the specific values of a 
and (5. Since \im.f l ^ f ± 00 g/ r (h) = ±00 and £/ r (0) = —3/2 < , it has always at least two roots, one positive 
and one negative. For the same reason, it cannot have two positive and two negative roots, since at each 
simple root the quartic function changes sign. Since the function srf (h) changes smoothly with a and /3, 
the boundaries between regions where there are five roots and regions where there are three roots are 
found enforcing that £/(h) has a multiple root. This condition is satisfied only by the points belonging 
to the parabolas (3 = c + a 2 and (3 = c_ a 2 , where c + = 1/4 and c_ is the only real root of the equation 
8 + 48 y — 435 y 2 + 676 y 3 = 0. The regions above the positive parabola and below the negative one have 
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Figure 7: phase space diagram for the number of asymptotic solutions 



only three roots, which are simple roots, while the regions between the two parabolas (except f3 = 0) have 
five roots, which are again simple roots. On the boundaries (3 = c±a 2 between the three-roots regions and 
the five-roots regions there are four roots, one of which is a root of multiplicity two. This is summarised 
in figure 7. 

We name the roots in the following way: the y = root is denoted as L. For the phase space points 
where there are just three roots, the positive root is denoted as C+ and the negative one as C- . For 
points in the five-roots regions, we adopt the following convention. Be («5, f3§) a point where there are five 
roots. In the same quadrant of the phase space, take another point (0:3,^3) where there are three roots, 
and a path ^€ which connects the two points. Following the path ^ ', two of the four non-zero roots of 
(0:5, f3§) smoothly flow to the non-zero roots of (03, ^3), and are denoted as C+ and C_ themselves. The 
other two non-zero roots of (0:5, /%), instead, disappear when (following ^) the boundary of the five-roots 
region is crossed, and are denoted as Pi and P2. We adopt the convention that |Pi| < |!F*2 1 - The definition 
is independent of the particular choice of the point (0:3, ^3) and of the path ^ used. A careful study of 
the asymptotic equation and of its derivatives permits to show that we have C < C + < Pi < P2 for 
a > and P2 < Pi < C_ < C + for a < 0. On the boundaries /3 = c± a 2 we have Pi = P2 = P. 

Inner solutions 

Suppose now that a solution of the quintic equation exists in a neighborhood of p = + (possibly not 
defined in p = 0) , and that it has a well defined limit when p — > + . Then such a solution cannot tend to 
zero as p — > + . In fact, for p 7^ we can divide (A.l) by A, and indicating x = p/p v we obtain 

x 3 Q^ 2 h 5 - (a 2 + 2/3) h 3 + 3ah 2 - ^h^j + 3(3 h 2 -1 = (A.6) 

If we had h — > as p — > + , then we would get —1 = when taking the p — > + limit in the previous 
equation. Therefore, lim p _ 5 . + h(p) / . In a neighbourhood of p = + , we can further divide the equation 
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above by /i 5 , obtaining the following equation in v = 1/h 

t,5 + ! X 3 „ 4 - 3 (0 + a x 3 ) v 3 + (a 2 + 2/3) x 3 t> 2 - ^(3 2 x 3 = (A.7) 

The permitted limiting values for v are then the roots of the equation obtained from the previous one 
setting x = 0, namely 

v 5 -3(3v 3 = (A. 8) 

For (3 > there are three roots, namely vq = 0, v + = +y / 3~/3 and v_ = — y/3~J3 ; for f3 < 0, instead, there 
is only the root v = 0. Therefore, the permitted limiting behaviours for h when p — > + are 

\h(p)\ -»• +oo (A.9) 

for /3 / 0, and 
only for /? > 0. 

A. 2 Leading behaviours 

Note that so far we have not proved that inner and asymptotic solutions exist, but just found the values 
that have to be the limit of these solutions if they exist. The existence and uniqueness of solutions can be 
proved applying the implicit function theorem (known also as Dini's theorem, see for example [19]) to the 
equation (A.l) around ^4 = and to the equation (A.7) around x = 0. One can then prove that, to each of 
the roots vq, v± of (A. 8), it is possible to associate a local solution of (3.3) defined in a neighbourhood of 
p = + , which are respectively the diverging inner solution D and the finite inner solutions F±. Likewise, 
it can be shown that it is possible to associate a local solution of (3.3) defined in a neighbourhood of 
p = +oo, to each of the simple roots of (A. 4). These solutions are the asymptotic solutions L, C+, C_, 
Pi and P2. In the case Pi = P2 = P, a separate analysis is needed. It can be shown that for a > and 
(3 = c + a 2 there are no local solutions of (3.3) which tend to P when p — > +00, and the same holds for 
a < and j3 = c_ a 2 . On the other hand, for a > and j3 = c_ a 2 there are two different local solutions 
of (3.3) which tend to P when p — > +00, and the same holds for a < and /3 = c + a 2 . Despite having the 
same limit for p — > +00, these two local solutions are different when A / 0: we then call Pi the solution 
which in absolute value is smaller, and P2 the solution which in absolute value is bigger. Therefore, 
on the boundaries between the three-roots-at-infinity regions and the five-roots-at-infinity regions, for 
a ^ 0, f3 = c± a 2 there are three asymptotic solutions of (3.3), while for a ^ 0, j3 = c T a 2 there are five 
asymptotic solutions of (3.3). 

We now turn to the discussion of the leading behaviours of the inner and asymptotic solutions. For 
the finite inner solutions F± and finite non-zero asymptotic solutions C± and Pi,2, the behaviour is 

h{p) = C + R{p) (A.ll) 

where C 7^ is their limiting value, and R is respectively such that lim p _ 5> o+ R = (inner solutions) and 
lim p ^ +00 R = (asymptotic solutions). For the asymptotic decaying solution L and the inner diverging 
solution D, a more detailed study is worthwhile. 

Asymptotic decaying solution L 
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Let's consider the solution L, which satisfies lim p ^ +00 h(p) = 0. Dividing the quintic equation (A.l) 
by h, we get 



lp 2 h 4 -(a 2 + 2p)h 2 + 3ah-l=(^) (j.-30h) (A.12) 

The left hand side has a finite limit when p — > +00, so the same has to hold for the right hand side: 
taking this limit in the equation above gives 

lim f^Yl = -| (A.13) 

p^+oo \ p h 2 K ' 



Hp) = -Uj) 3 + R(p) (A.14) 



which implies that 

with lmip^+oo p 3 R(p) = 0. 
Inner diverging solution D 

Let's consider now the solution D, which satisfies lim p ^ +\h(p)\ = +00. Dividing the equation (A. 7) 
by v 3 , one finds that 



v 



2 3/3= (-) 3 l(-^ 4 + 3Q " 3 -( a2 + 2 ^ 2 + ^ 2 ) (A.15) 



One more time, the left hand side has a finite limit when p — > + , so the same should hold for the right 
hand side. Therefore, the p — > + limit in the equation above gives 

lim (£-) 3 ± = -l (A.16) 

and so 

«(P) = -^|- + R(P) (A.17) 
V 2 Pv 

with lim p _ >0 + R(p)/p = 0. To understand the behaviour of the gravitational potentials (3.1)-(3.2) in this 
case, it is useful to calculate the next to leading order behaviour. In fact, it turns out that, after going back 
to h = 1/v, the leading behaviour precisely cancels the Schwarzschild-like contribution, so to understand 
if the gravitational potentials are finite at the origin it is essential to know how R behaves for very small 
radii. Inserting (A.17) into (A. 7) and dividing by x 5 , one obtains taking the limit p — > + that 

lim J^ = J-f« 2 + -/3 > \ (A.18) 

where x = p/p v . We have then 



\ 2 p v \p v J 



where 



N=±{* + §*) (A.20) 
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and lim p ^ +(TZ(p)/p 3 ) = 0. Finally, going back to the function h we get 



where 



V P P Pv 



M = l\lw( a2 + 3 2^ (A - 22) 



and ]ha p _ i . Q +(^(p) / p) = 0. It can be shown that in the special case a 2 + 3/3/2 = 0, the next to leading 
order term scales as p 2 instead of p, and that \i~m p ^Q+(&(p) / p 2 ) = 0. 

Therefore, we can conclude that in general the diverging inner solution D is such that 



h(p) = -^^+R(p) (A.23) 



where lim /9 _ 5>0 + (R(p) I P) ls finite. 
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